INSTABILITY OF THE LINEARIZED PROBLEM IN MULTIWAVE 
TOMOGRAPHY OF RECOVERY BOTH THE SOURCE AND THE SPEED 



PLAMEN STEFANOV AND GUNTHER UHLMANN 

Abstract. In this paper we consider the hnearized problem of recovering both the sound speed 
and the thermal absorption arising in thermoacoustic and photoacoustic tomography. We show 
that the problem is unstable in any scale of Sobolev spaces. 



1. Introduction 

In multiwave tomography, one sends one type of wave to the body of a patient, most often elec- 
tromagnetic (thermoacoustic tomography) or optical radiation (photoacoustic tomography) which 
interacts with the tissue, and measures the acoustic signal on the boundary generated by this 
interaction. This combines the high contrast of the incoming waves with the high resolution of 
the measured ultrasound ones. The mathematical model of the emitted ultrasound wave is the 
following. Let u solve the problem 

r {d^-c^A)u = in(0,r)xR", 

(1) < u\t=o = f, 

[ dtu\t=o = 0, 

where the sound speed c = c{x) > and T > are fixed. Assume that / and c — 1 are supported 
in O, where C R" is some bounded domain with a smooth convex boundary. The measurements 
are modeled by the operator 

(2) Ai/ := u|[o,T]x9Q- 

The first step in multiwave imaging is to recover / given Ai/. The speed c is usually assumed to 
be known but in practice, it is not. Then the natural question is whether we can recover both c 
and /. The answer is still unknown. A discussion of this problem, together with a partial local 
result stating that c can be determined up to a constant scaling can be found in |10] . David Finch 
observed a link between this problem and the transmission eigenvalues. 

The problem of recovery of /, given c has received a lot of attention in the past years. We refer 
to [IliaEliaElElElElElIIOlIiaillllSlIITlIISlI^ for some works in this direction. If c 

is constant, then the inversion of Ai is actually an integral geometry problem as well. For variable 
c, there is always uniqueness if T ^ 1, and there is stability if c is non-trapping and T 3> 1. We 
refer to [23l [20j for details. 

When an inaccurate speed is used for reconstruction by time reversal, the resulting images are 
distorted, and full of "artifacts", see for example the images in |12j or [10]. The mathematical 
structure of the artifacts is easy to understand: the forward map Ai is a Fourier Integral Operator 
(FIG) with a canonical relation given by the graph of the map determined by the geodesies rays 
from the interior to the boundary, see [23] for more details. The time reversal is another FIG with a 
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canonical relation given by the graph of the inverse of the former one. When the speed is different, 
we get that time reversal with a wrong speed is an FIO with a canonical relation not the diagonal 
(the graph of the identity) but the graph of the composition of the forward one and a backward 
one with a different speed. In particular, since the canonical relation of Ai and its reversal have 
two disconnected components, one can see double images of the same singularity, well visible in the 
examples in [T2], for instance. An algorithm to tune in the speed by maximizing the sharpness of 
the reconstructed / is proposed in [26j. This is related to the FIO description of the artifacts but 
good mathematical understanding of this algorithm is lacking. 

One of the ways to recover the speed is to take additional measurements and to recover c from 
travel times. This is the method proposed in [27J in thermoacoustic tomography. The travel time 
problem is stable under some geometric assumptions on c (c~^dx^ being a simple metric in the 
domain is enough) which are satisfied when c is close enough to a constant, in particular, see, e.g., 
[19] and the references there. Then / can be recovered stably as well. Additional data can be 
provided either by placing ultrasound sources around the body, or by placing passive absorbing 
(tissue imitating) objects around the body which become ultrasound sources by the thermo- or 
the photo-acoustic effect. We refer to the recent paper [12] for references and numerical and 
experimental implementations of that method. As can be expected, the results are very good. 
Simultaneous reconstruction of / and c aside from the above mentioned work [26j , have been tried 
with various success in [HU [301 [32] , for example. 

In this paper, we study the linearization 5Ai and we show that the latter is unstable. In partic- 
ular, we prove the following. 

Theorem 1. There is no stability estimate of the type 



si >0,S2> 0, regardless of si, S2. 

We also show that a conditional type of stability estimate cannot hold either, see Remark [TJ This 
suggests instability of the non-linear problem as well but does not imply it directly. Stability of the 
linearization in some Sobolev norms, even if not in the sharp ones, does imply (conditional Holder) 
stability of the non-linear problem [22]. The converse however is a much more delicate question. 
To prove the instability of the linearization, we show that the latter is a smoothing operator for 
{5f^5(?) belonging to an explicitly defined infinite dimensional linear space, see (|39p. We refer to 
section [H for more details. 

2. Preliminaries 

Notice first that c^A is formally self-adjoint w.r.t. the measure c~^dx. Given a domain U , and 
a function u{t,x), define the energy 



In particular, we define the space H£,{U) to be the completion of Cq'{U) under the Dirichlet norm 



It is easy to see that H£){U) C H^{U), if U is bounded with smooth boundary, therefore, H£,{U) 
is topologically equivalent to Hq{U). If [/ = R", this is true for n > 3 only. By the finite speed 



¥f\\Hni^n) + ¥c^\\H^riK) < C\\5Ai{5f,5c''}\ 



H=2 ' 
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of propagation, the solution with compactly supported Cauchy data always stays in even when 
n = 2. The energy norm for the Cauchy data [/i, /2], that we denote by || ■ ||^ is then defined by 

ll[/,/2]|||^= / (|V./lP + C-2|/2|2)dx. 

Ju 

This defines the energy space 

n{U)=HD{U)(BL\U). 
Here and below, = c^^dx). Note also that 

(4) \\f\\j,^ = {-c'Af,f)L.. 
The wave equation then can be written down as the system 

(5) ut = Pu, P = , P := c^A, 

where u = [u, ut] belongs to the energy space Ti. The operator P then extends naturally to a skew- 
selfadjoint operator on if c € and € In this paper, we will deal with either U = R" 
or [/ = Q. In the latter case, the definition of Hd(U) reflects Dirichlet boundary conditions. 

We generalize next the results in ^23j to the inverse problem with general Cauchy data (/i,/2) 
in ([T]) with g not necessarily zero. What we really need later is Proposition [2] only. Let u solve the 
problem 

r {df-c'^A)u = in(0,T)xR", 

(6) < u\t=o = fl, 

[ dtu\t=o = /2, 

where T > is fixed. Set f = [/i, /s]. Then for f G "H, we have u G C([0, T]; U). 

Assume that / is supported in 0, where Vt C R" is some smooth bounded domain. Set 

(7) Af := u|[o,T]x9C- 

The trace Af is well defined in C[q-^{^[Q,T\; H^^'^{dQ)) , where the subscript (0) indicates that we 
take the subspace of functions h so that h = for t = 0. For a discussion of other mapping 
properties, we refer to [H]. When f is restricted to functions in V., supported in a fixed compact 
/C C rj, then A is an FIO with a canonical relation of graph type, and maps f continuously into 
^(0) ([0'^] ^ see [23]. 



Given h, let v solve 



(8) 



{di - c^A)v = in (0, T) x n, 

'"\lo,T]xdn = h, 

v\t=T = 4>, 

dtv\t=T = 0, 



where cj) solves the elliptic boundary value problem 

(9) A0 = O, 4>\dn = h{T,-). 
Then we define the following pseudo-inverse 

(10) A/i := [w(0,-),t't(0, •)] inf^. 
By [11], 

A : /7fo)([0,r] -xdn)^H = H^{n) X L'^{n) 
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is a continuous map. Note that the mapping properties above allow us to apply A to Af only when 
f is compactly supported in but the theorem above shows that AA extends continuously to the 
whole %. 

Let T(rj) be the length of the longest geodesic in il, when c~'^dx^) is non-trapping. 

Theorem 2. Let (VL^c^'^dx^) he non-trapping, and let T > T{Q,). Then AA = Id—K., where K 
is compact in 7i{il.), and |[K||-^(q) < 1. In particular, /d — K is invertihle on H(f]), and A has an 
explicit left inverse of the form 

oo 

(11) f = K™A/i, h := Af. 

m=0 

Proof Let f G x C7~(0) first. Let w solve 

{df - c^A)w = in (0, T) x n, 

w\[o,T]xdn = 0, 



(12 



W\t=T = U\t=T - 
Wt\t=T = Ut\t=T, 



where u solves ([6]) with a given f € 7^. Let v be the solution of ([8]) with h = Af. Then v + w 
solves the same initial boundary value problem in [0, T] x 0, that u does (with initial conditions at 
t = T), therefore u = v + w. Restrict this to t = to get 

f = AAf + vi^(0, •)• 

Set 

Kf = w(0,-) = [wiO,-),wtiO,-)]. 

We will show now that K extends to a compact operator. Since T > T{Q), all singularities starting 
from leave ft at t = T. Therefore, u{T,-) and ut{T,-), restricted to fl, are C°°. Moreover, 
considered as linear operators of f , they are operators with smooth Schwartz kernels. Then so is 
(p, see dHl), by elliptic regularity. Therefore, the map 

n{Q)3i ^ [u{T,-)-^,ut{T,-)]en{n), 

defined a priori on Cq°{^}) x Cq^{^}) functions, extends to a compact one because it is an operator 
with smooth kernel on 0. Since the solution operator of (jl2p from t = T to t = is unitary in 7^(fi), 
we get that the map 7i{Q) 9 f i— >• ['w{0, •),wt{0, •)] G H£){i^)) is compact, too, as a composition of 
a compact and a bounded one. 

We know remove the smoothness restriction on f , and let it be any element in T-L. In what follows, 
(•, ■)hd{^) is inner product in i?£)(J7), see ([3]), applied to functions that belong to H^{i^) but 
maybe not to H^i^) (because they may not vanish on dQ). Set := u{T, ■). By ^ and the fact 
that = (/) on do., we get 

{u^ - (/>,(/>) = 0. 

Then 



Therefore, the energy of the initial conditions in ()12p satisfies the inequality 
(13) En{w, T) = \\u^ - m^^^^ + ||uf < En{u, T). 

Since the Dirichlet boundary condition is energy preserving, we get 

En{w,0) = En{w,T) < En{u,T) < Enn{u,T) = En{u,0) = ||f|||^. 
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Therefore, 

(14) \\Kf\\l^^^^=En{w,0)<\\{\\l^^^y 
We show next that actually the inequality above is strict, i.e., 

(15) ||Kf||^(f,) < i|f|l^(f,), f/0. 

Assume the opposite. Then for some f 7^ 0, all inequalities leading to (fill) are equalities. In 
particular, E^{w,T) = Ejin{u,T). Then 

n(T, x) = 0, for x ^ $7. 

By the finite domain of dependence then 

(16) n(t,x)=0 when dist(x,0) > |r - t|. 
One the other hand, we also have 

(17) u(t,x) = when dist(x,$7) > \t\. 
Therefore, 

(18) u{t, x) = when dist(x, dQ) > T/2, -T/2 < t < 3T/2. 

In [23] , we used the fact that u extends to an even function of t that is still a solution of the wave 
equation because /2 = there. Then one gets that (jlSp actually holds for \t\ < 3T/2. Then one 
concludes by Tataru's unique continuation theorem [25] that u = on [0,T] x f), therefore, / = 0. 
We also noted there that the time interval [0,T] is actually larger (twice as large) than what we 
need for the Neumann series to converge, see |24t l2Uj . where T > T{Q)/2 only. 

In the case under consideration, /2 does not necessarily vanish. We modify the arguments as 
follows. Prom John's theorem (equivalent to Tataru's |23l Theorem 2] in the Euclidean setting), 
we get that n = on [0, T] x \ Q. Then [221 Theorem 2] implies that n = for t = T/2 and all 
X. By energy preservation, f = 0. 

Now, one has 

(19) ||Kf||^(f,) < yA7||f||«(^,), f/0, 

where Ai is the largest eigenvalue of K*K. Then Ai < 1 by (fT5]) . □ 
Denote by 

B := (Id - K)-^A 
the left inverse of A constructed in Theorem [2l 

Proposition 1. 

B : i7j,)([0,r] xdn)^'H = H^{n) X L'^{n) 

is a continuous map. 

Proof. Note first, that A has the mapping properties above, by the results in |16J. Next, (Id — K)^^ 
is a bounded map in 7i by Theorem [2l □ 

Let Bi^2 be the components of B (that sends scalar functions to vector functions), i.e., B/i = 
{Bih, B2h). Let A12 be the components of A (that sends vector functions to scalar functions, 
i.e., Af = Ai/i + A2/2. We can think of B as a 2x1 matrix, and of A as an 1x2 matrix. Then 
BA[/i,/2] = [/i,/2], therefore, 

(20) BiAi = Id, B2A2 = Id, B1A2 = 0. B2A1 = 0. 

Set d^^g = g{s) ds and let Ad he the Dirichlet realization of A in U. We have the following. 
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Proposition 2. 



(21) 



dtki on HD{^)nH^ 

Ai, A2 = di^Ai onHoin), 
Ai(c2Ad)"\ onL^in), 



Proof. Let u solve ([6]) with f = [/i,0] G -D(P), see ([5]). Then solves the wave equation, too, 
with Cauchy data [0,c'^Afi] G n. This proves the first relation in (j2ip . Similarly, let u solve ([6]) 
with f = [0, /a] G -D(P). Then dtu solves the wave equation, too, with Cauchy data [/2,0]. This 
proves the second relation. 

For / G H]:){Q), we have S^Aa/ = Ai/ by what we just proved. Moreover, A2/ = for t = 0. 
This proves the third relation in (j2ip . 

Finally, let / G L^{n). Then dtAiic^Ao)'^ f = A2/ by the first identity in (HH. Moreover, 
Ai{c^ Ad)~^ f = for t = 0. This proves the fourth identity. □ 

Relations (f2T|) can be used to show that A2 has a left inverse based on the result in [23] only. 
Indeed, set B2 := Bidt- Then ^2^2 = BidtA2 = BiAi = Id on Hni^)- Analyzing the mapping 
properties of B'2 and A2 as in [23], one gets that the composition .BiAi is well defined and is a 
bounded operator on L^(il); that therefore has to be identity. 

3. Recovery of the speed c when / is known; The linearization SAi/dc^ w.r.t. the 



The non-linear problem of recovery of c when / is known, and its linearization were studied in 
detail by the authors in |21j . We showed there that if 6c^ is a priori supported in a compact subset 
/C, then the linearization SAi/dc^ is Fredholm. We also gave conditions for uniqueness for both the 
linear and the non-linear problem. The geometric requirement is that there exists a foliation of Cl 
by strictly convex surfaces. If n = 2, non-trapping implies that condition. In all dimensions, if c is 
close enough to a constant, that condition holds. In particular, if 



then the parts of the spheres |x| = R, R > 0, intersecting Cl form a foliation of surfaces convex w.r.t. 
the metric c~'^dx'^, see \20\ 121]. There is an another condition as well: we require that Af{x) ^ 0, 
Vx G K. This condition fits well into our analysis; it is an if and only if condition for the operator 
Qn below to be elliptic. 



We derive a representation of the difference Ai/— Ai/ first. If one is interested in the linearization 
only, the computations are not really simpler — one just has to drop the tilde in most of the 
formulas. 



SPEED 



(22) 



X ■ dc < c in Q 



4. Analysis of the linearized operator 6A 



1 



(23) 



Let (c,/), C 




Then 



A1/-A1/ 



(n — u) 



[0,T]xdQ' 



We have 



(24) 
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where w solves 



(25) 
Set 



- c'^A)w = (c2 - c2)An in (0,r) x R", 
w\t=o = 0, 
dtw\t=o = 0. 

P:=c'A, h = c-\??-^). 

ft 



and let Pd be the Dirichlet realization of P in U. Recall the notation df- = Jq g{s) ds. Denote 
by Ru with u = [ui, U2] the restriction of ui to [0, T] x dQ. We have 



T]xdQ = ^l U{s)[0,hPu{t-s,-)]ds. 



wl 

Take the t derivative to get 

(26) wt\^^^^^^g^ = A2hPf + W, W ■.= Rj^U{s)[0,hPut{t-s,-)]ds. 
Differentiate W to get 

(27) = R [ U{s) [0, hd^Pu{t - s, •)] ds, 

<Ji Jo 

where we used the fact that dtu = for t = 0, therefore dtAu = for t = 0. Differentiate one more 
time to get 

(28) —^W = K2hP'^f + R U{s)[i:),hdlPu{t- s,-)]ds, 

Ot Jq 

Relations ([26]) and 1^ show that W = 0{t'^) at t = 0. Therefore, 

(29) W = d^'^A2hP^f + d^'^R [ U{s)[0,hdfPu{t-s,-)]ds:=h+l2. 

Jo 

By Proposition [21 for the first term on the right we get 

h = AsP^i (hP^f 



The second term on the r.h.s. of (j29|) can be written in the form 

l2:=d^^R [ U{t- s)[0,hP^ut{s,-)]ds. 
Jo 

We claim that I2 = I21 where 

= /" U{t- s)[0, P^^hP^ut{s, •)] ds. 
Jo 

Indeed, I( = for t = 0, and direct differentiation shows that 

dtl'2 = R I U{t- s)[Pi^^hP'^ut{s, ■),{)] ds. 
Jo 

because RPj^^ = 0. Therefore, dtl'i = for t = as well. Differentiate one more time to get 
I2 = -^2- Therefore, 

W = K2Pj^^hP'^f + R [ U{s)[0,P^^hP^ut{t- s,-)]ds. 

Jo 
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Compare with ([26]) and repeat these arguments to get for any = 1, 2, . . . 

W =y^K2PD^hP^^^~f + R / U{s)[{),P^^+^hP^ut{t- s,-)]ds. 



We want to emphasize that in all those formulas, h is considered as an operator of multiplication, 

^^^hP^+^f = K2{PD' 



i.e., K2PD^hP^+^f = K2{PD^{hP^^^f)), etc. Combine this with ([MD, to get the following. 



(30) 5jAi/-Ai/ 



= k2P{f-f)+Y,K2PD^hP^+^f + R / U{s)%P^''+^hP''ut{t-sr)]ds. 

Since Aif — Aif = for t = 0, integrating w.r.t. t, and applying Proposition [21 we get 
(31) AiZ-Ai/ 

= AiCf - f) + y^KiP^^hP^~f + d;^R / t/(5)[0,P^^/iP^+int(t-s,-)]ds. 

Finally, using the arguments above, we write the last term on the right as 

R f U{s)[P^^-^hP^+^ut{t - s, •), 0] ds. 
Jo 

We therefore proved the following. 

Theorem 3. Let (c, /) and (c, /) be two pairs in C"^^ (n) x H'^'^+'^{n) n H^{n) with c> 0, c > 
in Q. Then 



N 



Ai/ - Ai/ =Ai(/ - /) + Ai ^ (c^Az,)-' (c"2(c2 _ c') {c'A)'f 



k=l 



(32) 



+ R 



U{s) [{c'An)-''-' {c~\c'-c'){c'Af^%{t-s,.)) ,0 



ds. 



We define the linearization 5Ai{5/, (5c^} at (c, /) as the derivative at e = of Ai with speed 
= c + e5c and source fe = f + s^f- 

Corollary 1. For any N = 1,2, . . . , the linearized operator 5Ai{5f, Sc^} at (c, /) has the form 

(33) dAi{6f,dc^} = Ai{6f + Qnc~^6c^) +Rnc'^Sc^, 
where 

N 

(34) QNh = (c^Ab)-' (hic'Aff] 



k=l 



and 
(35) 



RNh = rJ^U{s) [(c^Ad)"^"^ (/i(c2A)^+^nt(t-s,-)) ,0 



ds. 
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Proof of TheoremUl Note that for (c, /) G C°°, the operator i?jv is smoothing 2N + 2 degrees, i.e., 

(36) Rn : H^{1C) ^ H'+^^+\[0, T] x On), 

and Qn is a ^'DO (in the interior of Q) with principal symbol A/(a;)|^|~^. So for supported in 
a compact set in 17, we can write 

(37) SAi{6f, 6c^} = Ai {6f + (c^^ A/II^I^^ ^ ^^2^ ^ iioo(^c2, 

where i?oo is smoothing, and l.o.t. stands for "lower order terms", i.e., for a pseudo-differential 
operator of order —3, which we can always assume to have a proper support. 

This result is not surprising. All singularities of the kernel of 6Ai are of conormal type, at 
dist(a;, y) = t, where dist is the distance in the metric c~^dx^. This can be seen from the represen- 
tation 

(38) dAi{6c^,6f} = Ai6f + R [ U{s)[0, Au{t - s, ■)6c'^]ds. 

Jo 

The operator Qtv can be explained by those singularities. Next, Rj^ depends on the smooth part 
of the kernel, and in particular, the behavior of the kernel inside that geodesic ball. 

An important observation is that if the expression in the parentheses in (133p vanishes (and this 
is an explicit condition on 6f and (5c^), then the linearization is very smooth. If for a moment we 
ignore Roo, we get a kernel of infinite dimension. Indeed, given Sc^, compactly supported in Q, 
we can always find 6f so that 6f + (c~^A/|D|~^ -|- l.o.t.) 6c^ = 0. If we are given a compactly 
supported (in Q) function 6f, and if A/ on supp(5/, we can find 6c^ at least away from a 
finitely dimensional space, so that the expression above vanishes as well. The effect of -Roo will be 
the following. Even though we still do not know whether we can determine both 6f and 6c^, we 
can claim that even of we could, the problem would be unstable. 

More precisely, given > 1, let Vat be the linear space 

(39) Vn = {{5f,h) e Hoin) X L\lCy, 5f + QNh = 0}. 
Then for the linearization 6Ai{6f,6c'^} we have 

(40) Vn 3 {5f,c-H^) ^ 5Ai{5f,5^] G hI^^+\%T] x dn). 

This can easily be generalized to {5h, f) belonging to negative Sobolev spaces. We then complete 
the proof of Theorem [1] by a well known argument, see, e.g., [22j or the remark below. □ 

Remark 1 . It is easy to see also that 6Ai does not satisfy the following conditional stability estimate 
either: 

(41) \\Sf\\Hnin) + \\Sc^\\H^Hic)<C\\6Ai{6f,5c'}\\l,^, < ^ < 1, 
under the condition 

(42) \\Sf\\H^s(n) + \\Sc'\\H^si,c)<A 

regardless of how we choose si, S2, S3 and ^ > 0. Indeed, fix (/> G C^, supported in the interior of 
IC. Set 

where A > 0, S3 > and uj G S"""^ are fixed. Let 6 fx be such that {5f\,h\) G Vn, A^ ^ 1. Then 
(I42p is satisfied. On the other hand, with = {1 + hx)(? , 

||Mi{5/A,5ci}||^,, = \\RNhx\\H^2 < CnX''-''-^^-^ 
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see (|36|) . while 

Therefore we get a contradiction with (j4ip for ^ 1. 
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